# Figure_A09.R

# Part of the replication archive for 
#
#   Bullock, John G. 2020. "Education and Attitudes toward Redistribution in
#   the United States." British Journal of Political Science 50.


# This file produces Figure A9 in the appendix to the article: "Associations 
# of outcomes with age, period, and cohort."


library(Bullock, lib.loc = c(.libPaths(), 'packageLibrary'))       # for qw()
library(car)           # for Recode()
library(dplyr)         # for %>%, group_by(), left_join(), mutate()
library(lattice)
library(grid)
library(gridExtra)     # for grid.arrange()
library(RColorBrewer)
library(tibble)        # for rownames_to_columns()
library(tools)         # for toTitleCase()

source("IV_setup.R")

COMPLEX_PANEL_NAMES   <- qw('cohortAge cohortPeriod')
COMPLEX_PANEL_OUTCOME <- 'eqwlth'
COMPLEX_PANEL_SPAN    <- 2/3     # Default is 2/3. Higher values will produce more smoothing. 
namesDF <- data.frame(
  varNames = qw('eqwlth goveqinc guarantee.7pt govt.health.7pt helppoor welfare'),
  dfNames  = qw('GSS.df GSS.df   ANES.df       ANES.df         GSS.df   GSS.df'),
  varName_pretty = c(
    'redist. (1)', 'redist. (2)', 'guarantee', 'health', 'help poor', 'welfare'),
  stringsAsFactors = FALSE)



##############################################################################
# PRELIMINARIES FOR PRINTING THE FIGURE 
##############################################################################
dirOutput        <- 'float_output/'  # where the figures go
filenameStem     <- 'Figure_A09' 
PDFtitle         <- 'Figure A9: Associations between age, period, cohort, and the outcomes of interest'
PS_width         <- 9.5
PS_height        <- 12.0
postscriptBackground <- 'transparent'
plotLineColor    <- brewer.pal(6, 'Dark2')  
  plotLineColor[1:2] <- plotLineColor[2:1]
plotLineWidth    <- 1
panelBorderCol   <- 'black'
panelBorderWidth <- .3
panelNumberCol   <- 'gray37' 
panelLayout        <- c(3, 1)  # columns, then rows
panelLayoutComplex <- c(2, 1)
xBetween         <- 2.00       # space between columns
xBetweenComplex  <- 6.00       # space between columns
yBetween         <- 2.30       # space between rows
yBetweenComplex  <- 4.30       # space between rows
baseCexSize      <- 1          
xAxisTextSize    <-  .8*.87    # cex
yAxisTextSize    <-  .8*.87    # cex
xLabSize         <- baseCexSize * .85
yLabSize         <- xLabSize
axisTickSize     <- .4
yAxis            <- list(
  draw        = FALSE,
  labels      = qw('.4 .5 .6 .7'), 
  at          = c(.4, .5, .6, .7),
  limits      = c(.375, .725),
  tck         = c(axisTickSize, 0),
  col         = panelBorderCol,
  cex         = yAxisTextSize,
  alternating = TRUE,
  relation    = 'same',  
  rot         = 0,
  axs         = 'i')
yAxisComplex     <- list(
  draw        = FALSE,
  labels      = c('.4', '', '.45', '', '.5'), 
  at          = seq(.4, .5, by = .025),
  limits      = c(.375, .525),
  tck         = c(axisTickSize, 0),
  col         = panelBorderCol,
  cex         = yAxisTextSize,
  alternating = TRUE,
  relation    = 'same',  
  rot         = 0,
  axs         = 'i')

# SET PANEL WIDTH AND HEIGHT
panelWidth         <- list(1.50, "inches")         
panelHeight        <- list(2.25, "inches")  
panelWidthComplex  <- list(2.25, "inches")          
panelHeightComplex <- list(3.08, "inches")  



##############################################################################
# ADD DECADE-GROUPED AGE AND COHORT VARIABLES TO ANES.DF AND GSS.DF
##############################################################################
ANES.df$ageGrouped <- cut(
  x              = ANES.df$age,
  breaks         = c(-100, seq(20, 80, by=10), 150),
  right          = FALSE,
  ordered_result = TRUE)
GSS.df$ageGrouped <- cut(
  x              = GSS.df$age,
  breaks         = c(-100, seq(20, 80, by=10), 150),
  right          = FALSE,
  ordered_result = TRUE)
ANES.df$yearIntGrouped <- cut(
  x              = ANES.df$yearInt,
  breaks         = c(-100, seq(1979, 1999, by=10), 2099),
  dig.lab        = 4,
  ordered_result = TRUE)
GSS.df$yearIntGrouped <- cut(
  x              = GSS.df$yearInt,
  breaks         = c(-100, seq(1979, 1999, by=10), 2099),
  dig.lab        = 4,
  ordered_result = TRUE)
ANES.df$yearYoungGrouped <- cut(
  x              = ANES.df$yearYoung,
  breaks         = c(-100, seq(1919, 1999, by=10), 2099),
  dig.lab        = 4,
  ordered_result = TRUE)
GSS.df$yearYoungGrouped <- cut(
  x              = GSS.df$yearYoung,
  breaks         = c(-100, seq(1919, 1999, by=10), 2099),
  dig.lab        = 4,
  ordered_result = TRUE)
ANES.df$yearYoungGrouped20 <- cut(
  x              = ANES.df$yearYoung,
  breaks         = c(-100, seq(1919, 1999, by=20), 2099),
  dig.lab        = 4,
  ordered_result = TRUE)
GSS.df$yearYoungGrouped20 <- cut(
  x              = GSS.df$yearYoung,
  breaks         = c(-100, seq(1919, 1999, by=20), 2099),
  dig.lab        = 4,
  ordered_result = TRUE)



##############################################################################
# CREATE DATA FRAME FOR PLOTTING OF SIMPLE PANELS
##############################################################################
# We start by creating the data frame for the first three panels, which just
# show bivariate relations between the six outcomes and either age, period, or
# cohort.  [2018 08 01]
#
# Each buildDF() function builds a data frame containing all of the data that 
# we need for one outcome variable. We use as.data.frame() because ANES.df is 
# actually a tibble, which causes problems with the code that I use to extract 
# the att and educ vectors.  
buildDF_ungroupedData <- function (varName, ...) {
  dfName    <- namesDF$dfNames[namesDF$varNames == varName] 
  attDF     <- as.data.frame(get(dfName))
  att       <- attDF[, varName]

  buildDF_ungroupedData_helper <- function (varName, att, xAxisVar, predictorValueName) {
    data.frame(
      varName          = varName, 
      meanOutcome      = att, 
      xAxisVar         = xAxisVar,
      xAxisValue       = attDF[, predictorValueName],
      stringsAsFactors = FALSE)
  }
  age    <- buildDF_ungroupedData_helper(varName = varName, att = att, 'age',    'age')
  period <- buildDF_ungroupedData_helper(varName = varName, att = att, 'period', 'yearInt')
  cohort <- buildDF_ungroupedData_helper(varName = varName, att = att, 'cohort', 'yearYoung')
    
  bind_rows(age, period, cohort) %>%
    .[complete.cases(.), ] %>%
    mutate(xAxisVar = factor(xAxisVar, levels = qw('age period cohort')))
}

APC_descriptives.df <- mapply(
    FUN = buildDF_ungroupedData,
    namesDF$varNames,
    MoreArgs = list(grouping = 'NONE'),
    SIMPLIFY = FALSE) %>%
  Reduce(bind_rows, .)  


# ADD INTEGER VALUES FOR PREDICTORS
# The predictors -- i.e., the x-axis values, age or period or cohort -- are
# specified in the xAxisValue vector. If it is a character vector, it will be 
# converted to a factor variable by xyplot(). If the data are not grouped, 
# there will be no problem because xAxisValue will be a numeric vector. But 
# if the data are grouped (by year or by 10-year period), the conversion 
# to a factor will screw things up. For example, by creating a single factor, 
# xyplot() will be putting "[1910, 1920)" and "1930" on the same scale. And 
# a given year value, like 1970, will have a single level in the factor, even
# though it should have one level when period is on the x-axis, and quite 
# another when cohort is on the x-axis. These problems will lead to screwy 
# plotting.
#   We solve the problem by making each predictor value correspond to an 
# integer. For example, when the predictor is age, predictor value "(-100, 20]" 
# corresponds to 1. But when the predictor is period, "1970" corresponds to 1.
#   Creation of integers for the period (interview year) variable is tricky.
# 1970 is the first year in which any question was asked, so it gets a 1. But
# in 1971, no outcome question was asked. Even so, 1971 should be counted as 2,
# and 1972 (the next year in which outcome questions were asked) should be 
# counted as 3. If I simply count 1972 as 2 -- skipping "unused" years -- the
# panels that plot period will be misleading.
makePeriodIntegers <- function(myDF) { 
  myDF %>% 
    ungroup %>%
    filter(., xAxisVar == 'period') %>%  
    select(., xAxisValue) %>%
    unlist %>%
    as.integer %>%
    Bullock::rescale(., c(1, length(min(.):max(.))))
}

APC_descriptives.df$xAxisValueInt <- APC_descriptives.df$xAxisValue  


# ADD LINE LABELS TO DATA FRAME
# The labels are the "varName_pretty" variable in namesDF.
APC_descriptives.df <- left_join(
  x = APC_descriptives.df,
  y = select(namesDF, -dfNames),
  by = c("varName" = "varNames")) 


# ORDER COLUMNS
# This re-ordering of columns doesn't affect plotting. It's for the  
# convenience of those who will read this file. It also makes the ordering of 
# columns in APC_descriptives.df correspond to that of 
# APC_descriptives_complex.df.
APC_descriptives.df <- APC_descriptives.df %>%
  select(
    xAxisVar, xAxisValue, xAxisValueInt, 
    varName, varName_pretty, meanOutcome, everything())



##############################################################################
# CREATE DATA FRAME FOR PLOTTING OF COMPLEX PANELS
##############################################################################
# Whereas the "simple" panels plot only bivariate relationships, the "complex" 
# panels do more. Each complex panel plots data for a single outcome. One of 
# { age, period, cohort } is along the x-axis. And the different lines in the 
# panel correspond to different values of another { age, period, cohort } 
# variable. For example, we might plot age along the x-axis, and have different
# lines within the panel correspond to different cohort groupings. In all 
# cases, the y-axis indicates the value of the outcome, just as with the simple
# panels.  
buildDF_complex <- function (varName, panelNames, ...) {
  dfName <- namesDF$dfNames[namesDF$varNames == varName] 
  attDF  <- as.data.frame(get(dfName))
  att    <- attDF[, varName]
  
  buildDF_complex_helper <- function (
    varName, 
    att, 
    attDF, 
    xAxisVarname, 
    xAxisVarnamePretty, 
    curveVarname,
    curveVarnamePretty) {    
      data.frame(
        meanOutcome      = att,
        xAxisValue       = attDF[, xAxisVarname],
        curveVarValue    = as.character(attDF[, curveVarname]),
        varName          = varName, 
        xAxisVar         = xAxisVarnamePretty, 
        curveVar         = curveVarnamePretty,
        panel            = paste0(xAxisVarnamePretty, toTitleCase(curveVarnamePretty)),
        stringsAsFactors = FALSE) %>%
      add_column(curveVarValuePretty = Recode(
        var     = .$curveVarValue, 
        recodes = '"(-100,1919]"="1910-19"; "(-100,1979]"="1970-79"; "(1919,1929]"="1920-29"; "(1929,1939]"="1930-39"; "(1939,1949]"="1940-49"; "(1949,1959]"="1950-59"; "(1959,1969]"="1960-69"; "(1969,1979]"="1970-79"; "(1979,1989]"="1980-89"; "(1989,1999]"="1990-99"; "(1999,2099]"="2000-08"; "(-100,20]"="18-20"; "[20,30)"="20-29"; "(20,30]"="21-30"; "[30,40)"="30-39"; "(30,40]"="31-40"; "[40,50)"="40-49"; "(40,50]"="41-50"; "[50,60)"="50-59"; "(50,60]"="51-60"; "(60,70]"="61-70"; "(70,80]"="71-80"; "(80,150]"="81+"')) %>%
      select(xAxisVar, xAxisValue, curveVar, curveVarValue, curveVarValuePretty, varName, meanOutcome, panel)       
    }
  
  agePeriodMeans    <- buildDF_complex_helper(varName, att, attDF, 'age',       'age',    'yearIntGrouped',   'period')
  ageCohortMeans    <- buildDF_complex_helper(varName, att, attDF, 'age',       'age',    'yearYoungGrouped', 'cohort')  
  periodCohortMeans <- buildDF_complex_helper(varName, att, attDF, 'yearInt',   'period', 'yearYoungGrouped', 'cohort')  
  cohortAgeMeans    <- buildDF_complex_helper(varName, att, attDF, 'yearYoung', 'cohort', 'ageGrouped',       'age')
  cohortPeriodMeans <- buildDF_complex_helper(varName, att, attDF, 'yearYoung', 'cohort', 'yearIntGrouped',   'period')

  # CHOOSE THE PANELS THAT WILL BE PLOTTED IN THE FIGURE
  # We use panelNamesFactor both to set the order of the panels (via 
  # "mutate(panel = ...") and to filter out data on panels that we won't use
  # (via "mget(.)"). 
  panelNamesFactor  <- factor(panelNames, levels = panelNames)

  # CREATE THE DATA FRAME
  paste0(panelNamesFactor, 'Means') %>%
    mget(., inherits = TRUE) %>%
    bind_rows %>%
    .[complete.cases(.), ] %>%
    mutate(panel = factor(panel, levels = levels(panelNamesFactor)))
}

APC_descriptives_complex.df <- buildDF_complex(
  varName    = COMPLEX_PANEL_OUTCOME, 
  panelNames = COMPLEX_PANEL_NAMES)


# FILTER OUT SOME DATA
# The panels aren't intelligible when they contain 10 or 15 lines -- for, say,
# 10 or 15 cohorts. I could collapse the cohorts into 20-year periods. Or I 
# can limit the display to only those cohorts that have the most data. With
# this code, I take the latter approach. I am trying to limit each panel to 
# five lines.
#
# Note: loess() can't plot the eqwlth~period relationship for those in the 
# (1999, 2099] cohort, because there are only five observations for this 
# cohort: not five individual cases, but five periods in which the question 
# was asked. This is a further reason to remove this particular cohort from 
# the data, at least from the periodCohort panel.  
cohortsToRemove   <- qw('(-100,1919] (1919,1929] (1929,1939] (1939,1949] (1989,1999] (1999,2099]')
ageGroupsToRemove <- qw('[-100,20) [60,70) [70,80) [80,150)')
APC_descriptives_complex.df <- APC_descriptives_complex.df %>%
  filter(!(panel%in%qw('ageCohort periodCohort') & curveVarValue%in%cohortsToRemove)) %>%
  filter(!(panel%in%qw('cohortAge periodAge')    & curveVarValue%in%ageGroupsToRemove))  


# TURN curveVarValuePretty INTO A FACTOR WITH ORDERED LEVELS
# Lattice is going to turn it into a factor later if we don't do it now. And
# we need to get the levels ordered correctly -- all age groups together, and
# all period groups together -- if the colors are to come out right. 
# [2018 08 25]
if (all(COMPLEX_PANEL_NAMES == c('cohortAge', 'cohortPeriod'))) {
  APC_descriptives_complex.df$curveVarValuePretty <- factor(
    x      = APC_descriptives_complex.df$curveVarValuePretty,
    levels = qw('1970-79 1980-89 1990-99 2000-08 20-29 30-39 40-49 50-59'))
} 



##############################################################################
# SET UP X AXES
##############################################################################
xAxis_labels <- list(
  qw('25 35 45 55 65 75'),    
  qw('1976 1986 1996 2006'),
  qw('1925 1950 1975 2000'))
xAxis_at <- lapply(xAxis_labels, as.integer)

xAxis <- list(
  draw        = TRUE, 
  labels      = xAxis_labels,
  at          = xAxis_at,
  tck         = c(axisTickSize, 0), 
  col         = panelBorderCol, 
  cex         = xAxisTextSize,
  alternating = 1, 
  relation    = 'free', 
  axs         = 'i')  # axs='i' means that there is no padding around the xlimits

# X-AXIS LIMITS
xAxis$limits <- list(
  c(18, 82),
  NULL,
  c(1920, 2005))


# X AXES FOR COMPLEX PANELS
xAxisComplex_labels <- list(
    agePeriod    = qw('25 35 45 55 65 75'),  
    ageCohort    = qw('25 35 45 55 65 75'),
    periodCohort = qw('1980 1990 2000 2010'),
    cohortAge    = qw('1925 1950 1975 2000'),
    cohortPeriod = qw('1925 1950 1975 2000')) %>%
  .[levels(APC_descriptives_complex.df$panel)]  # filter our list elements that apply to panels that we won't plot

xAxisComplex_at <- lapply(xAxisComplex_labels, as.integer)
xAxisComplex_limits <- list(
    agePeriod    = c(18, 82),          
    ageCohort    = c(18, 82),          
    periodCohort = NULL,            
    cohortAge    = c(1922, 2003),
    cohortPeriod = c(1922, 2003)) %>%
  .[levels(APC_descriptives_complex.df$panel)] # %>%  # filter our list elements that apply to panels that we won't plot

make_XAxisComplex <- function (limits, labels, at) {
  list(
    draw        = TRUE, 
    lim         = limits,
    labels      = labels,
    at          = at,
    tck         = c(axisTickSize, 0), 
    col         = panelBorderCol, 
    cex         = xAxisTextSize,
    alternating = 1, 
    relation    = 'free', 
    axs         = 'i')  
} 

# xAxisComplex is a traditional list. It's for a single figure that has 
# multiple panels.
#   xAxisComplex_indivPanels is a list of lists. Each "interior" list is a 
# traditional "xAxis" object, fit for a one-panel figure.  
xAxisComplex <- make_XAxisComplex(xAxisComplex_limits, xAxisComplex_labels, xAxisComplex_at)
xAxisComplex_indivPanels <- mapply(
  FUN      = make_XAxisComplex, 
  limits   = xAxisComplex_limits,
  labels   = xAxisComplex_labels, 
  at       = xAxisComplex_at, 
  SIMPLIFY = FALSE)



##############################################################################
# PANEL FUNCTION FOR SIMPLE PANELS
##############################################################################
APC_descriptives.panel <- function(x, y, ...) { 

  # PLOT THE DATA
  panel.xyplot(x, y, ...)

  
  # AXIS LABELS AND TICKS
  # I cannot manipulate clipping as easily as I want. And it is therefore 
  # a bit difficult to create axis labels and ticks that are outside the 
  # panels: I want clipping to be on so that data (i.e., plotted lines) don't
  # extend outside the panels, but I want to be able to draw labels and ticks
  # outside the panels. Judicious use of par.settings(clip) and 
  # trellis.par.set("clip", list(panel="off")) doesn't help here. 
  #   There are two solutions. The Procrustean solution is to filter the data
  # so that they correspond to the xlimits. The other solution, which I adopt 
  # here, is to draw internal ticks in this panel function, external ticks and
  # labels after printing the plot.
  panel.axis(
    side     = 'top', 
    at       = xAxis$at[[panel.number()]], 
    labels   = NULL,
    outside  = FALSE,
    half     = FALSE,
    tck      = axisTickSize)
  panel.axis(
    side     = 'left', 
    at       = yAxis$at, 
    labels   = NULL,
    outside  = ifelse(panel.number()==1, TRUE, FALSE),
    half     = FALSE,
    tck      = axisTickSize)
  panel.axis(
    side     = 'right', 
    at       = yAxis$at, 
    labels   = NULL,
    outside  = ifelse(panel.number()==3, TRUE, FALSE),
    half     = FALSE,
    tck      = axisTickSize)
  
  
  # LINE LABELS
  # all label coordinates provide upper left-hand corner of label
  simplePanel_1_lineLabels <- list(
    name  = sort(namesDF$varNames),  # order in which colors are chosen
    label = c('redist. (1)', 'redist. (2)', 'health\ncare', '',  'help\npoor', 'welfare'),
    x     = unit(c(26,       44,            65,             00,  55,           37),  'native'),
    y     = unit(c(.41,      .6175,         .4475,          100, .5425,        .69), 'native')) 
  simplePanel_2_lineLabels <- list(
    name  = namesDF$varNames,
    label = c('redist. (1)',  '',  'health care',  'guaranteed\nstd. of living',  'help\npoor',  'welfare'),
    x     = unit(c(1979.5,    00,  1972,           1972,                          55,            1990),  'native'),
    y     = unit(c(.430,      00,  .51,            .6125,                         .545,          .695),   'native')) 
  simplePanel_3_lineLabels <- list(
    name  = namesDF$varNames,
    label = c('redist. (1)',  'redist. (2)',  '',    '',    'help\npoor',  'welfare'),
    x     = unit(c(1925,      1944,           1972,  1972,  1925,          1943),  'native'),
    y     = unit(c(.430,      .62,            .52,   .6,    .53,           .690),  'native')) 
  simplePanel_labelLists <- list(simplePanel_1_lineLabels, simplePanel_2_lineLabels, simplePanel_3_lineLabels)
  grid.text(
    label = simplePanel_labelLists[[panel.number()]]$label,
    x     = simplePanel_labelLists[[panel.number()]]$x,
    y     = simplePanel_labelLists[[panel.number()]]$y,
    just = c('left', 'top'),
    gp    = gpar(cex=.75*yLabSize, col=plotLineColor, lineheight=.89))    
}



##############################################################################
# PANEL FUNCTION FOR COMPLEX PANELS
##############################################################################
APC_descriptives_complex.panel <- function(x, y, groups, panelName, span, ...) { 

  # PLOT THE DATA
  panel.superpose(x, y, groups, ...)
  
  # INTERIOR AXIS LABELS AND TICKS
  # Cannot plot exterior labels or ticks because I can't adjust clipping well
  # within the function. So I add the exterior labels and ticks after I print
  # the figure.
  #   This xAxisComplex$at[[panel.number()]] solution works only if I am  
  # plotting a single figure that has multiple panels. It won't work if I plot
  # multiple one-panel figures.
  panel.axis(
    side     = 'top', 
    at       = xAxisComplex$at[[panel.number()]], 
    labels   = NULL,
    outside  = FALSE,
    half     = FALSE,
    tck      = axisTickSize)
  panel.axis(
    side     = 'right', 
    at       = yAxisComplex$at, 
    labels   = NULL,
    outside  = FALSE,
    half     = FALSE,
    tck      = axisTickSize)

    
  # PANEL NUMBERS IN UPPER RIGHT-HAND CORNER OF EACH PANEL
  grid.polygon(
    x  = c(.875, .9975, .9975, .8750, .875), 
    y  = c(.900, .9000, .9975, .9975, .900), 
    gp = gpar(fill='gray87', col='gray87'))  
  grid.text(
    label = panel.number() + 3, 
    x     = .94, 
    y     = .950, 
    just  = 'center', 
    gp    = gpar(font=2, cex=1.15, col=panelNumberCol))

  
  # LINE LABELS 
  if (span == .8) {
    ageCohort_lineLabels <- list(
       label = c(paste0('1950\U00AD', '59'),   paste0('1960\U00AD', '69'),  paste0('1970\U00AD', '79'),  paste0('1980\U00AD', '89')),
       x     = unit(c(52.5,                    65,                          53,                          45.00),   'native'),
       y     = unit(c(.51,                     .4650,                       .4275,                       .46275),  'native')) 
    periodCohort_lineLabels <- list(
      label = c(paste0('1950\U00AD', '59'),   paste0('1960\U00AD', '69'),  paste0('1970\U00AD', '79'),  paste0('1980\U00AD', '89')),
      x     = unit(c(1997.75,                 2004,                        1984.50,                     1987.375),  'native'),
      y     = unit(c(.51,                     .4825,                       .435,                        .39),       'native')) 
    cohortAge_lineLabels <- list(
      label = c(paste0('20\U00AD', '29'),  paste0('30\U00AD', '39'),  paste0('40\U00AD', '49'),  paste0('50\U00AD', '59')),
      x     = unit(c(1987,                 1990,                      1980,                      1926.85), 'native'),
      y     = unit(c(.4025,                .4335,                     .4625,                     .44),     'native'))
    cohortPeriod_lineLabels <- list(
      label = c(paste0('1970\U00AD', '79'),   paste0('1980\U00AD', '89'),  paste0('1990\U00AD', '99'),  paste0('2000\U00AD', '08')),
      x     = unit(c(1929,                    1981,                        1925.75,                     1946),  'native'),
      y     = unit(c(.400,                   .382,                        .4815,                        .510),  'native'))
  } else {
    ageCohort_lineLabels <- list(
       label = c(paste0('1950\U00AD', '59'),   paste0('1960\U00AD', '69'),  paste0('1970\U00AD', '79'),  paste0('1980\U00AD', '89')),
       x     = unit(c(52.5,                    65,                          53,                          45.00),   'native'),
       y     = unit(c(.51,                     .4650,                       .4275,                       .46275),  'native')) 
    periodCohort_lineLabels <- list(
        label = c(paste0('1950\U00AD', '59'),   paste0('1960\U00AD', '69'),  paste0('1970\U00AD', '79'),  paste0('1980\U00AD', '89')),
        x     = unit(c(1997.75,                 2004,                        1984.50,                     1987.375),  'native'),
        y     = unit(c(.51,                     .4825,                       .435,                        .39),       'native')) 
    cohortAge_lineLabels <- list(
      label = c(paste0('20\U00AD', '29'),  paste0('30\U00AD', '39'),  paste0('40\U00AD', '49'),  paste0('50\U00AD', '59')),
      x     = unit(c(1967,                 1990,                      1980,                      1926.85), 'native'),
      y     = unit(c(.415,                .4300,                     .468,                     .4375),     'native'))
    cohortPeriod_lineLabels <- list(
      label = c(paste0('1970\U00AD', '79'),   paste0('1980\U00AD', '89'),  paste0('1990\U00AD', '99'),  paste0('2000\U00AD', '08')),
      x     = unit(c(1929,                    1984.5,                        1925.25,                     1944),  'native'),
      y     = unit(c(.400,                   .382,                        .4805,                        .515),  'native'))    
  }    
  lineLabel_list <- list(
    ageCohort    = ageCohort_lineLabels, 
    periodCohort = periodCohort_lineLabels,
    cohortAge    = cohortAge_lineLabels,
    cohortPeriod = cohortPeriod_lineLabels)
  
  currentPanelName <- panelName[panel.number()]  
  grid.text(
    label = lineLabel_list[[currentPanelName]]$label,
    x     = lineLabel_list[[currentPanelName]]$x,
    y     = lineLabel_list[[currentPanelName]]$y,
    just = c('left', 'top'),
    gp    = gpar(cex=.75*yLabSize, col=plotLineColor, lineheight=.89))      
}
  
  

##############################################################################
# CREATE THE PLOTS
##############################################################################
# When only one line is plotted per panel, plot.line governs the appearance of 
# the line.  When more than one line is plotted per panel, superpose.line 
# governs the appearance of the lines.
trellis.par.set(
  axis.components = list(     # distance between axis ticks and tick labels
    left  = list(pad1 = .50),
    right = list(pad1 = .50)),  
  axis.line = list(
    alpha = 1, 
    col   = panelBorderCol, # to eliminate panel border, set col to "white"
    lty   = 1, 
    lwd   = panelBorderWidth),
  superpose.line = list(col = plotLineColor)
)
noVerticalPadding_theme <-
 list(
   layout.heights = list(
    top.padding       = 0,
    main.key.padding  = 0,
    key.axis.padding  = 0,
    axis.xlab.padding = 0,
    xlab.key.padding  = 0,
    key.sub.padding   = 0,
    bottom.padding    = 0),
  layout.widths = list(
    left.padding      = 0,
    key.ylab.padding  = 0,
    ylab.axis.padding = 0,
    axis.key.padding  = 0,
    right.padding     = 0))


# CREATE PLOT OF SIMPLE PANELS
# These are the panels that simply show bivariate relations: attitude vs. 
# one of { age, period, cohort }. There are six lines in each panel: one for 
# each attitude.
APC_descriptives_plot <- xyplot(
  meanOutcome ~ xAxisValueInt | xAxisVar,
  groups   = varName,
  data     = APC_descriptives.df, 
  type     = 'smooth',  # 'smooth', 'l', 'p' for points, 'b' or 'o' for both line and points
  panel    = APC_descriptives.panel,
  as.table = TRUE,      # Plot panels top to bottom. Affects only multi-row figures.
  layout   = panelLayout,
  between  = list(y=yBetween, x=xBetween),                   
  strip    = FALSE,
  scales   = list(x=xAxis, y=yAxis),
  xlab     = '',
  ylab     = '',
  par.settings = list(
    clip = list(panel = "on"))
)


# CREATE PLOT OF COMPLEX PANELS
# makeComplexPanel() can be used to plot all panels in a single figure, or to 
# plot each panel as its own figure.  
makeComplexPanel <- function(complexDF_name, panelName, xAxisComplex) {
  complexDF <- get(complexDF_name)
  xyplot(
    meanOutcome ~ xAxisValue | panel,
    groups    = curveVarValuePretty,
    data      = filter(complexDF, panel %in% panelName),
    panelName = panelName,
    type      = 'smooth',  # 'smooth', 'l', 'p' for points, 'b' or 'o' for both line and points
    span      = COMPLEX_PANEL_SPAN,  # default is 2/3
    strip     = FALSE,
    panel     = APC_descriptives_complex.panel,
    layout    = panelLayoutComplex, 
    as.table  = TRUE,
    between  = list(y=yBetweenComplex, x=xBetweenComplex),                   
    scales    = list(x=xAxisComplex, y=yAxisComplex),
    xlab      = '',
    ylab      = '',
  )
}
APC_descriptives_complex_plot <- makeComplexPanel(
  complexDF_name = 'APC_descriptives_complex.df',
  panelName      = COMPLEX_PANEL_NAMES,         # e.g., "agePeriod"
  xAxisComplex   = xAxisComplex)
APC_descriptives_complexIndiv_plots <- mapply(  # returns a list of one-panel plots
  FUN            = makeComplexPanel, 
  complexDF_name = 'APC_descriptives_complex.df',
  panelName      = COMPLEX_PANEL_NAMES,
  xAxisComplex   = xAxisComplex_indivPanels,
  SIMPLIFY       = FALSE)



##############################################################################
# PRINT AND AUGMENT THE PLOTS
##############################################################################
# PRINT PLOT OF SIMPLE PANELS
print(
  x            = APC_descriptives_plot, 
  panel.width  = panelWidth, 
  panel.height = panelHeight)  


# ADD EXTERIOR LABELS AND TICKS TO SIMPLE PANELS, THEN CREATE GROB
xAxisLabels_simple <- c(
  'age', 
  'year of interview\n(\U201Cperiod\U201D)', 
  paste0('year turned 14\n(\U201C', 'cohort\U201D)'))  
for (panelNum in 1:3) {
  trellis.focus("panel", panelNum, 1, clip.off = TRUE, highlight = FALSE)
  panel.axis(
    side     = ifelse(panelNum==1, 'left', 'right'), 
    at       = yAxis$at, 
    labels   = if (panelNum==2) NULL else yAxis$labels,
    outside  = ifelse(panelNum==2, FALSE, TRUE),
    half     = FALSE,
    tck      = axisTickSize,
    text.cex = yAxisTextSize)
  grid.text(
    label = xAxisLabels_simple[panel.number()], 
    x     =  .5, 
    y     = -.133,
    just  = qw('center top'),
    gp    = gpar(font=1, cex=xLabSize, col=panelBorderCol, lineheight=.925), 
    rot   = 0)
  if (panelNum==1) {
    grid.text(
      label = 'average attitude', 
      x     = -.1933, 
      y     = .5, 
      gp    = gpar(font=1, cex=xLabSize, col=panelBorderCol), 
      rot   = 90)
  }
  
  # We draw the panel numbers (in the upper right-hand corner of each panel) 
  # here, rather than in the panel function, because we want to ensure that the 
  # ticks are behind the panel numbers rather than in front of them.
  grid.polygon(
    x  = c(.8575, .9975, .9975, .8575, .8575), 
    y  = c(.8925, .8925, .9975, .9975, .8925), 
    gp = gpar(fill='gray87', col='gray87', lineend='square', linejoin='mitre'))  
  grid.text(
    label = panel.number(), 
    x     = .93, 
    y     = .945, 
    just  = 'center', 
    gp    = gpar(font=2, cex=.9375, col=panelNumberCol))

  # We re-draw the panel border so that the polygon for the panel number 
  # doesn't overlap the border in any way.  [2018 08 20]
  grid.rect(gp = gpar(fill='transparent', lwd=panelBorderWidth))
  trellis.unfocus()
}
APC_descriptives_grob <- grid.grab(wrap = TRUE)


# PRINT PLOT OF COMPLEX PANELS
# We start by limiting plotLineColor to the number of curves (i.e., lattice
# groups) that we plot in each panel. Doing so makes it easy to match the line
# colors (selected by lattice from plotLineColor) and the label colors 
# (selected by me in the panel function). Without this command, values from 
# plotLineColor would be recycled in weird ways, leading to mixed-up colors.
trellis.par.set(
  superpose.line = list(col = plotLineColor[1:4]))

print(
  x            = APC_descriptives_complex_plot,
  panel.width  = panelWidthComplex, 
  panel.height = panelHeightComplex)  


# ADD EXTERIOR LABELS AND TICKS TO COMPLEX PANELS, THEN CREATE GROB
xAxisNames <- if (COMPLEX_PANEL_NAMES[2] == 'periodCohort') { 
  c('age', 'year of interview (\U201Cperiod\U201D)')
} else if (COMPLEX_PANEL_NAMES[2] == 'cohortPeriod') {
  paste0('year turned 14\n(\U201C', 'cohort\U201D)') %>%
    rep(., 2)
} else {
  sub('[A-Z].*', '', COMPLEX_PANEL_NAMES)
}

for (panelNum in 1:2) {
  trellis.focus("panel", panelNum, 1, clip.off = TRUE, highlight = TRUE)
  panel.axis(
    side     = 'left', 
    at       = yAxisComplex$at, 
    labels   = yAxisComplex$labels,
    outside  = TRUE,
    half     = FALSE,
    tck      = axisTickSize,
    text.cex = xAxisTextSize)
  grid.text(
    label = xAxisNames[panelNum], 
    x     =  .5, 
    y     = -.09,
    just  = qw('center top'),
    gp    = gpar(font=1, cex=xLabSize, col=panelBorderCol, lineheight=.925), 
    rot   = 0)
   grid.text(
     label = 'average attitude\n(only for first redistribution item)', 
     x     = -.19,
     # label = 'average attitude',
     # x     = -.17,
     y     = .5, 
     gp    = gpar(font=1, cex=xLabSize, col=panelBorderCol, lineheight=.9), 
     rot   = 90)
  trellis.unfocus()
}
APC_descriptives_complex_grob <- grid.grab()



##############################################################################
# ARRANGE PLOTS AND PRINT THE FIGURE  
##############################################################################
pdf(
  file   = paste0(dirOutput, filenameStem, '.pdf'), 
  width  = PS_width, 
  height = PS_height, 
  paper  = "special", 
  title  = PDFtitle,
  bg     = postscriptBackground)
  # cbPalette <- qw('#999999 #E69F00 #56B4E9 #009E73 #F0E442 #0072B2 #D55E00 #CC79A7')  # color-blind palette


# CREATE GRID LAYOUT
rectTransparent <- rectGrob(.5, .5, gp = gpar(col = 'transparent', fill = 'transparent'))
grid.arrange(
  grobs = list(APC_descriptives_grob, rectTransparent, APC_descriptives_complex_grob),
  ncol = 1,
  widths = NULL,
  heights = unit(c(panelHeight[[1]] + .5, .675, panelHeightComplex[[1]] + .5), 'in'))

dev.off()


# CROP THE PDF AND CLEAN UP INTERMEDIATE FILES
if (!(Sys.which('pdfcrop')=='' | Sys.which('perl')=='')  |  'pdfcrop' %in% installed.packages()[, 'Package']) {  # if "pdfcrop" is installed 
  system(paste(
    if (Sys.info()['sysname']=='Windows') paste(Sys.getenv('Comspec'), '/c ') else '',
    'pdfcrop', 
    paste0(dirOutput, filenameStem, '.pdf'), 
    paste0(dirOutput, filenameStem, '_crop.pdf')))
  file.remove(paste0(dirOutput, filenameStem, '.pdf'))
  if (interactive()) { 
    shell.exec(paste0(normalizePath(dirOutput), filenameStem, '_crop.pdf'))
  }
}